#include "CHBDYGAREACalculator.h"
#include "Utility/Algebra/Tools.h"

#include "DataStructure/Input/Element/CQUAD4Data.h"
#include "DataStructure/Input/Property/PCONVData.h"
#include "DataStructure/Input/Solver/ThermodynamicsConstant.h"

#include "DBManager/IService/Input/Node/INodeService.h"
#include "DBManager/IService/Input/Subcase/ISubcaseService.h"
#include "DBManager/IService/Input/Material/IThermalMaterialService.h"
#include "DBManager/IService/Input/Property/IPCONVService.h"

#include "DBManager/IService/Input/Element/IElementService.h"
#include "DBManager/IService/Input/Subcase/ISubcaseService.h"
#include "DBManager/IService/Input/Material/IThermalMaterialService.h"
#include "DBManager/IService/Input/Property/IPCONVService.h"
#include "DBManager/IService/Input/Load/ICONVService.h"
#include "DBManager/IService/Input/Constraint/ITEMPService.h"
#include "DBManager/IService/Input/Property/IPHBDYService.h"
#include "DBManager/IService/Input/Load/IQBDY3Service.h"

#include "SGFem/Integration/IntegrationFormula/Gauss2D2X2P.h"
#include "SGFem/Pre/Assembly/Assembly.h"
#include "SGFem/Pre/Assembly/ModelInfoQuery.h"

using SG::Algebra::Point;

using SG::Algebra::Vector3D;
using SG::Algebra::Matrix;
using SG::Algebra::identity;
using SG::Algebra::vstack;
using SG::Algebra::det3X3;
using SG::Algebra::calcTriArea;

using SG::DataStructure::Common::Status;
using SG::DataStructure::Common::Id;
using SG::DataStructure::Common::Real;
using SG::DataStructure::FEM::ThermalMaterial;
using namespace SG::DataStructure::FEM;
using namespace SG::DBManager;
using SG::FEMSolver::AssemblyEigen;

using SG::FEMSolver::Element::Gauss2D2X2P;
using SG::FEMSolver::Element::CHBDYGAREACalculator;
std::size_t CHBDYGAREACalculator::s_nNode4 = 4;
std::size_t CHBDYGAREACalculator::s_nNode3 = 3;
std::size_t CHBDYGAREACalculator::s_nDim = 3;

void CHBDYGAREACalculator::Initialize (SG::DataStructure::FEM::ElementTypes type)
{
    m_type = type;
    const auto& integ = m_integForm.GetPoints ();
    m_shpFunc.Compute (integ);
}

/// @brief 返回边界单元节点数量
void CHBDYGAREACalculator::NumofNodes(const Id eleType, _OUT std::size_t& nNode) const
{
    if (eleType == 5)
    {
        nNode = s_nNode4;
    }
    else if (eleType== 4)
    {
        nNode = s_nNode3;
    }
    else
    {
        // 丢出错误并终止程序
        throw std::runtime_error("Error: The boundary area element should be AREA3 or AREA.");
    }
}

/// @brief 使用求逆犯规法计算积分点Jacobian,并且组装应变矩阵`B`
void CHBDYGAREACalculator::ComputeStrainMatrix(std::size_t iGauss, const Matrix& coord, _OUT Matrix& B, _OUT Real& j) const
{
    const auto& DNDxi = m_shpFunc.GetShapeFunctionDerivative(iGauss);
    auto JT = DNDxi * coord;
    j = JT(0, 0) * JT(1, 1) - JT(0, 1) * JT(1, 0);
    auto JTinv = SG::Algebra::inverse2X2 (JT);
    B = JTinv * DNDxi;
}

/// @brief 计算积分点Jacobian
void CHBDYGAREACalculator::ComputedetJ(std::size_t iGauss, const Matrix& coord, _OUT Real& j) const
{
    const auto& DNDxi = m_shpFunc.GetShapeFunctionDerivative(iGauss);
    auto JT = DNDxi * coord;
    Matrix one(1, 3, {1.0, 1.0, 1.0});
    Matrix tmp = vstack(one, JT);
    j = abs(det3X3(tmp));
}

void CHBDYGAREACalculator::ComputeConvStiffness (const Matrix& coords, const ThermalMaterial& material, _OUT Matrix& ke2)
{
    const Real h =  material.m_H;
    
    if (coords.GetRowSize() == s_nNode3)
    {
        Real area;
        Vector3D normal;
        area = calcTriArea(coords, normal);

        ke2 = Matrix (3, 3, {2.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0});
        ke2 = ke2 * (h * area / 12.0);
    }
    else if (coords.GetRowSize() == s_nNode4)
    {
        Real j;
        const auto nGauss = m_integForm.GetPointSize();     
        for (std::size_t iGauss = 0; iGauss < nGauss; ++iGauss)
        {
            ComputedetJ (iGauss, coords, j);
            const Real tmp = j * m_integForm.GetPoint(iGauss).m_weight;
            auto& DN = m_shpFunc.GetShapeFunction(iGauss);
            auto DNT = DN.Transpose ();
            ke2  += DN * DNT * tmp;  
        }

        ke2 = ke2 * h; 
    }
}

void CHBDYGAREACalculator::ComputeConvStiffness
    (
    const SG::Algebra::Matrix& coords,
    const std::vector<Real>& h,
    _OUT SG::Algebra::Matrix& ke2
    )
{
    if (h.size() != 4 && h.size() != 3) {
        // 丢出错误并终止程序
        throw std::runtime_error("Error: The radiating boundary temperature point of the area element should be 3 or 4.");
    }
    
    if (coords.GetRowSize() == s_nNode3)
    {
        Real hsum =  (h[0] + h[1] + h[2]) / 3.0;
        Real area;
        Vector3D normal;
        area = calcTriArea(coords, normal);

        ke2 = Matrix (3, 3, {2.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0});
        ke2 = ke2 * (hsum * area / 12.0);
    }
    else if (coords.GetRowSize() == s_nNode4)
    {
        Real j;
        const auto nGauss = m_integForm.GetPointSize();     
        for (std::size_t iGauss = 0; iGauss < nGauss; ++iGauss)
        {
            
            ComputedetJ (iGauss, coords, j);
            const Real tmp = j * m_integForm.GetPoint(iGauss).m_weight;
            auto& DN = m_shpFunc.GetShapeFunction(iGauss);
            auto DNT = DN.Transpose ();
            Real hsum = h[0] * DN(0, 0) + h[1] * DN(1, 0) + h[2] * DN(2, 0) + h[3] * DN(3, 0);
            ke2  += DN * DNT * tmp * hsum;  
        }
    }
}

void CHBDYGAREACalculator::ComputeRadStiffness
    (
    const SG::Algebra::Matrix& coords, 
    const SG::DataStructure::Common::Real& Emissivity,
    const std::vector<SG::DataStructure::Common::Real>& T,
    const SG::DataStructure::Common::Real& Text,
    _OUT SG::Algebra::Matrix& ke4
    )
{
    if (T.size() != 4 && T.size() != 3) {
        // 丢出错误并终止程序
        throw std::runtime_error("Error: The radiating boundary temperature point of the area element should be 3 or 4.");
    }

    //hr最大长度为4，如果是三角形边界单元，则最后一位可以不管
    std::vector<Real> hr(s_nNode4);
    for (std::size_t i = 0; i < T.size(); ++i)
    {
        hr[i] = STEFAN_BOLTZMANN_NUMBER * Emissivity * (T[i] * T[i] + Text * Text) * (T[i] + Text);
    }

    ComputeConvStiffness(coords, hr, _OUT ke4);
}

void CHBDYGAREACalculator::ComputeConvLoad
    (
    const SG::Algebra::Matrix& coords, 
    const SG::DataStructure::FEM::ThermalMaterial& material,
    const SG::DataStructure::Common::Real& Text,
    _OUT SG::Algebra::Matrix& p3
    )
{
    Real h =  material.m_H;

    if (coords.GetRowSize() == s_nNode3)
    {
        Real area;
        Vector3D normal;
        area = calcTriArea(coords, normal);

        p3 = Matrix (3, 1, {1.0, 1.0, 1.0});
        p3 = p3 * (h * Text * area / 3.0);

    }
    else if (coords.GetRowSize() == s_nNode4)
    {
        Real j;
        const auto nGauss = m_integForm.GetPointSize();     
        for (std::size_t iGauss = 0; iGauss < nGauss; ++iGauss)
        {
            ComputedetJ (iGauss, coords, j);
            const Real tmp = j * m_integForm.GetPoint(iGauss).m_weight;
            auto& DN = m_shpFunc.GetShapeFunction(iGauss);
            p3  += DN * tmp;  
        }
        p3 = p3 * h * Text;
    }
}

void CHBDYGAREACalculator::ComputeConvLoad
    (
    const SG::Algebra::Matrix& coords, 
    const std::vector<Real>& h,
    const SG::DataStructure::Common::Real& Text,
    _OUT SG::Algebra::Matrix& p3
    )
{
    if (h.size() != 4 && h.size() != 3) {
        // 丢出错误并终止程序
        throw std::runtime_error("Error: The radiating boundary temperature point of the area element should be 3 or 4.");
    }

    if (coords.GetRowSize() == s_nNode3)
    {
        Real hsum =  (h[0] + h[1] + h[2]) / 3.0;
        Real area;
        Vector3D normal;
        area = calcTriArea(coords, normal);

        p3 = Matrix (3, 1, {1.0, 1.0, 1.0});
        p3 = p3 * (hsum * Text * area / 3.0);
    }
    else if (coords.GetRowSize() == s_nNode4)
    {
        Real j;
        const auto nGauss = m_integForm.GetPointSize();     
        for (std::size_t iGauss = 0; iGauss < nGauss; ++iGauss)
        {
            ComputedetJ (iGauss, coords, j);
            const Real tmp = j * m_integForm.GetPoint(iGauss).m_weight;
            auto& DN = m_shpFunc.GetShapeFunction(iGauss);
            Real hsum = h[0] * DN(0, 0) + h[1] * DN(1, 0) + h[2] * DN(2, 0) + h[3] * DN(3, 0);
            p3  += DN * tmp * hsum;  
        }
        p3 = p3 * Text;
    }
}

void CHBDYGAREACalculator::ComputeHeatFluxLoad
    (
    const SG::Algebra::Matrix& coords, 
    const SG::DataStructure::Common::Real& q0,
    _OUT SG::Algebra::Matrix& p2
    )
{
    // if (q0.size() != 4) {
    //     // 丢出错误并终止程序
    //     throw std::runtime_error("Error: The Heat Flux boundary points of the AREA element should be 4.");
    // }
    if (coords.GetRowSize() == s_nNode3)
    {
        Real area;
        Vector3D normal;
        area = calcTriArea(coords, normal);

        p2 = Matrix (3, 1, {1.0, 1.0, 1.0});
        p2 = p2 * (q0 * area / 3.0);
    }
    else if (coords.GetRowSize() == s_nNode4)
    {
        Real j;
        const auto nGauss = m_integForm.GetPointSize();     
        for (std::size_t iGauss = 0; iGauss < nGauss; ++iGauss)
        {
            ComputedetJ (iGauss, coords, j);
            const Real tmp = j * m_integForm.GetPoint(iGauss).m_weight;
            auto& DN = m_shpFunc.GetShapeFunction(iGauss);
            Real qsum = q0 * DN(0, 0) + q0 * DN(1, 0) + q0 * DN(2, 0) + q0 * DN(3, 0);
            p2  += DN * tmp * qsum;  
        }
    }
}

void CHBDYGAREACalculator::ComputeHeatFluxMatries 
    (
        const SG::DBManager::DBServiceFactorySP dbServiceFactory,
        const SG::DataStructure::FEM::CHBDYGData eleData,
        const SG::DataStructure::Common::Real& q0,
        const SG::DataStructure::Common::Id sid
    )
{
    std::size_t nNode = 0;
    NumofNodes(eleData.m_TYPE, _OUT nNode);

    auto& assembly = AssemblyEigen::GetInstance ();
    
    static const DofSet DOFs {DOFDirections::DOF_Ux};

    std::vector<std::size_t> nodeIndexes(nNode);
    std::vector<Id> dispCoordIds(nNode);
    std::vector<int> dofGlobalIndex (nNode);
    std::vector<Point> points(nNode);

    const auto& nodeService = dbServiceFactory -> get<INodeService>();

    // Node coordinate
    ModelInfoQuery::QueryNode (nodeService, eleData.m_g, nNode, _OUT points, _OUT nodeIndexes, _OUT dispCoordIds);
    
    
    Matrix coords(nNode, s_nDim);
    for (std::size_t iNode = 0; iNode < nNode; ++iNode)
    {
        coords(iNode, 0) = points[iNode].m_x;
        coords(iNode, 1) = points[iNode].m_y;
        coords(iNode, 2) = points[iNode].m_z;
    }
    assembly.GetGlobalDOFId (nodeIndexes, DOFs, _OUT dofGlobalIndex);

    //热流荷载贡献
    Matrix p2(nNode, 1);
    ComputeHeatFluxLoad(coords, q0, _OUT p2);
    assembly.AccumulateLoad(sid, dofGlobalIndex, p2);
}

void CHBDYGAREACalculator::ComputeConvMatries 
    (
        const SG::DBManager::DBServiceFactorySP dbServiceFactory,
        const SG::DataStructure::Common::Id sid
    )
{
    auto& assembly = AssemblyEigen::GetInstance ();
    
    static const DofSet DOFs {DOFDirections::DOF_Ux};

    ThermalMaterial material;
    CHBDYGData eleData;
    PCONVData property;
    CONVData convData;
    //TEMPData tempData;
    //QBDY3Data heatFluxData;

    //std::vector<QBDY3Data> heatFluxDataArray;
    std::vector<TEMPData> tempDataArray;
    
    const auto& chbdygService = dbServiceFactory ->get<IElementService>();
    const auto& nodeService = dbServiceFactory -> get<INodeService>();
    const auto& convService = dbServiceFactory ->get<ICONVService>();
    const auto convSize = convService-> Size();

    for (std::size_t index{0}; index < convSize; ++index)
    {
        /** 从数据库获取对应的对流数据 */
        convService->FindByIndex (index, _OUT convData);

        if (Status::STATUS_Failure == chbdygService->FindById (convData.m_id, _OUT eleData))
        {
            continue;
        }
        
        // 确定边界单元的节点数目
        std::size_t nNode = 0;
        NumofNodes(eleData.m_TYPE, _OUT nNode);

        std::vector<std::size_t> nodeIndexes(nNode);
        std::vector<Id> dispCoordIds(nNode);
        std::vector<int> dofGlobalIndex (nNode);
        std::vector<Point> points(nNode);

        // Node coordinate
        ModelInfoQuery::QueryNode (nodeService, eleData.m_g, nNode, _OUT points, _OUT nodeIndexes, _OUT dispCoordIds);
        
        Matrix coords(nNode, s_nDim);
        for (std::size_t iNode = 0; iNode < nNode; ++iNode)
        {
            coords(iNode, 0) = points[iNode].m_x;
            coords(iNode, 1) = points[iNode].m_y;
            coords(iNode, 2) = points[iNode].m_z;
        }
        assembly.GetGlobalDOFId (nodeIndexes, DOFs, _OUT dofGlobalIndex);

        //获取材料参数，这部分怎么写？
        dbServiceFactory ->get<IPCONVService>()->FindById (eleData.m_id, _OUT property);
        dbServiceFactory ->get<IThermalMaterialService>()->FindById (property.m_MID, _OUT material);

        //获取对流换热系数和热交换温度
        dbServiceFactory ->get<IPCONVService>()->FindById (convData.m_PCONID, _OUT property);
        dbServiceFactory ->get<IThermalMaterialService>()->FindById (property.m_MID, _OUT material);
        
        //对流刚度矩阵贡献
        Matrix ke2(nNode, nNode);
        ComputeConvStiffness(coords, material, _OUT ke2);
        assembly.AccumulateStiffness(dofGlobalIndex, ke2);

        //对流荷载贡献
        Matrix p3(nNode, 1);
        Real Text;
            
        dbServiceFactory ->get<ITEMPService>()->FindBySId(sid, _OUT tempDataArray);
        for (auto& tempData : tempDataArray)
        {
            if (tempData.m_G == convData.m_TA[0])
            {
                Text = tempData.m_T;
                break;
            }
        }
        ComputeConvLoad (coords, material, Text, _OUT p3);
        assembly.AccumulateBCLoad(sid, dofGlobalIndex, p3);
    }
}

void CHBDYGAREACalculator::ComputeElementMatries (const SG::DBManager::DBServiceFactorySP dbServiceFactory)
{
    auto& assembly = AssemblyEigen::GetInstance ();
    
    static const DofSet DOFs {DOFDirections::DOF_Ux};

    // SG::DataStructure::FEM::CQUAD4Data eleData;
    // PCONVData property;
    ThermalMaterial material;
    CHBDYGData eleData;
    PCONVData property;
    CONVData convData;
    TEMPData tempData;
    QBDY3Data heatFluxData;
    // ThermalMaterial material;
    // CONVData convData;
    std::vector<QBDY3Data> heatFluxDataArray;
    std::vector<TEMPData> tempDataArray;
    
    std::vector<std::size_t> nodeIndexes(s_nNode4);
    std::vector<Id> dispCoordIds(s_nNode4);
    std::vector<int> dofGlobalIndex (s_nNode4);
    std::vector<Point> points(s_nNode4);
    
    const auto& nodeService = dbServiceFactory -> get<INodeService>();
    const auto& chbdygService = dbServiceFactory ->get<IElementService>();
    const auto eleSize = chbdygService->GetElementSize (m_type);
    for (std::size_t index{0}; index < eleSize; ++index)
    {
        /** 从数据库获取对应的单元数据 */
        chbdygService->FindByIndex (index, _OUT eleData);

        // Node coordinate
        ModelInfoQuery::QueryNode (nodeService, eleData.m_g, s_nNode4, _OUT points, _OUT nodeIndexes, _OUT dispCoordIds);
        Matrix coords(s_nNode4, s_nDim);
        for (std::size_t iNode = 0; iNode < s_nNode4; ++iNode)
        {
            coords(iNode, 0) = points[iNode].m_x;
            coords(iNode, 1) = points[iNode].m_y;
            coords(iNode, 2) = points[iNode].m_z;
        }
        assembly.GetGlobalDOFId (nodeIndexes, DOFs, _OUT dofGlobalIndex);

        //获取材料参数，这部分怎么写？
        dbServiceFactory ->get<IPCONVService>()->FindById (eleData.m_id, _OUT property);
        dbServiceFactory ->get<IThermalMaterialService>()->FindById (property.m_MID, _OUT material);
        
        Matrix ke(s_nNode4, s_nNode4);
        ComputeConvStiffness(coords, material, _OUT ke);
        assembly.AccumulateStiffness(dofGlobalIndex, ke);


        Id sid = 1;   //这个地方后期后改一改
        //获取对流换热系数和热交换温度
        if (Status::STATUS_Success == dbServiceFactory ->get<ICONVService>()->FindById (eleData.m_id, _OUT convData))
        {
            // for (auto& convData : convDataArray)
            // {
            dbServiceFactory ->get<IPCONVService>()->FindById (convData.m_PCONID, _OUT property);
            dbServiceFactory ->get<IThermalMaterialService>()->FindById (property.m_MID, _OUT material);
        
            //对流刚度矩阵贡献
            Matrix ke2(s_nNode4, s_nNode4);
            ComputeConvStiffness(coords, material, _OUT ke2);
            assembly.AccumulateStiffness(dofGlobalIndex, ke2);

            //对流荷载贡献
            Matrix p3(s_nNode4, 1);
            Real Text;
            dbServiceFactory ->get<ITEMPService>()->FindBySId(sid, _OUT tempDataArray);
            for (auto& tempData : tempDataArray)
            {
                if (tempData.m_G == convData.m_TA[0])
                {
                    Text = tempData.m_T;
                    break;
                }
            }
            ComputeConvLoad (coords, material, Text, _OUT p3);
            assembly.AccumulateLoad(sid, dofGlobalIndex, p3);
        }

        // sid = 2;
        // //获取热流参数
        // if (Status::STATUS_Success == dbServiceFactory ->get<IQBDY3Service>()->FindBySId (sid, _OUT heatFluxDataArray))
        // {   
        //     for (auto& heatFluxData : heatFluxDataArray)
        //     {
        //         if (eleData.m_id == heatFluxData.m_EID)
        //         {
        //             //热流荷载贡献
        //             Matrix p2(s_nNode4, 1);
        //             Real q0 = heatFluxData.m_Q0;
        //             ComputeHeatFluxLoad(coords, q0, _OUT p2);
        //             assembly.AccumulateLoad(sid, dofGlobalIndex, p2);
        //         }
        //     }
        //     continue;
        // }

    }
}


// 后处理数据
void CHBDYGAREACalculator::DataRecovery (const SG::DBManager::DBServiceFactorySP dbServiceFactory, int subcaseId, SG::DataStructure::Common::Real time, const std::vector<SG::DataStructure::Common::Real>& globalDisplacement)
{
}

// 计算全局坐标系下单元应力刚度矩阵
void CHBDYGAREACalculator::ComputeStressStiffness (const SG::DBManager::DBServiceFactorySP dbServiceFactory, int staticSubId)
{
}
